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This article discusses the basic properties of the path integral 
method for continuum fermions, focusing on the restricted path in- 
tegral (RPIMC) approach. The RPIMC is a simulation technique 
for fermion systems that extends the previously developed bosonic 
methods to degenerate many- fermion systems, avoiding the fermion 
sign problem by restricting the paths. For general fermion systems, 
the restriction needed for exact results is not known, so ansatz for 
the nodes of the fermion density matrix are used. Even with this 
approximation, the finite temperature Monte Carlo method is ex- 
pected to be even more accurate and powerful than corresponding 
zero-temperature Monte Carlo methods. 
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I. INTRODUCTION 

The lectures at this school have concerned various 
types of simulations. We have heard about various at- 
tempts to perform quantum dynamical calculations, but 
what I want to tell you about are recent attempts to 
simulate the equilibrium properties of fermion systems. 
This remains one of the most challenging problems in the 
field and one of the most important both from a theoret- 
ical and practical point of view, given the pervasiveness 
of fermions in nature and intractability of research on 
many- fermion systems. 



Much progress has been made in recent years in 
performing simulations of quantum many-body boson 
and distinguishable-particle systems using imaginary- 
time path integrals. As first shown by Feynman (1953), 
and as we will later describe, the many-body density ma- 
trix can be obtained using classical-statistical methods 
on polymer-like systems. The density matrix of a many- 
body system at a temperature ksT = /3 _1 can be written 
as an integral over all paths {Rt}- 



p(R ,Rp;p) = ^Y l (±lf I dR t exp(-S[R 



])■ 



(1.1) 

The path R(t) begins at VRo and ends at Rp, and V is a 
permutation of particle labels. For N particles, the path 
is in 3N dimensional space: R t — (ri t ,r2t . . . rjyt). The 
upper sign is to be used for bosons and the lower sign for 
fermions. For nonrelativistic particles interacting with a 
potential V(R), the action of the path, S[Rt], is given 
by: 



f f 

S[Rt] = / 
Jo 



dt 



dR{t) 



hdt 



+ V{Rt) 



(1.2) 



Thermodynamic properties, such as the energy, are re- 
lated to the diagonal part of the density matrix, so that 
the path returns to its starting place or to a permutation 
V of its starting place after a "time" /?. 

Since the imaginary-time action S[R t ] is a real func- 
tion of the path, for boltzmannons and bosons the in- 
tegrand is nonnegative. It can thus be interpreted as 
a probability of an equivalent classical system and the 
action as the classical potential energy of a "polymer." 
To perform Monte Carlo calculations of the integrand, 
one makes imaginary time discrete, so that one has a fi- 
nite (and hopefully small) number of time slices and thus 
a classical system of N particles in M time slices; an 
equivalent NM particle classical system of "polymers." 
If the path integral is performed by a simulation method, 
such as a generalization of Metropolis Monte Carlo or 
with Molecular Dynamics, one can obtain essentially ex- 
act results for previously unsolved problems such as the 
properties of liquid 4 He at temperatures near the super- 
fluid phase transition, the exchange frequency in quan- 
tum crystals, and quantum particles immersed in classi- 
cal systems. (Ceperley, 1995) 
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Note that in addition to sampling the path, the permu- 
tation is also sampled. This is equivalent to allowing the 
ring polymers to connect in different ways. This macro- 
scopic "percolation" of the polymers is directly related 
to superfluidity as Feynman (1953) first showed. Any 
permutation can be broken into permutation cycles, i. 
e. into 2-, 3-, ... exchange cycles. Superfluid behavior 
can occur at low temperature when the probability of 
exchange cycles on the order of the system size is non- 
negligible. Figure 1 shows a trace of a typical path of 
six "superfluid" helium atoms in a 2D periodic square. 
I have recently written a detailed article describing the 
path integral theory of Bose superfluids and how one car- 
ries out the Monte Carlo calculations (Ceperley, 1995). 




FIG. 1. The extended trace of six 4 He atoms at a temperature 
of 0.75 K. The dashed square represents the periodic bound- 
ary conditions. The six filled circles represent the beginning 
of a path for each atom. Three of the atoms are involved 
in an exchange which winds around the boundary in the x 
direction. The paths have been Fourier smoothed for clarity. 



energy of the uniform electron gas using Diffusion Monte 
Carlo (Ceperley and Alder, 1980). In the fixed-node ap- 
proach, one solves the Schroedinger equation with the 
boundary condition that the many-body wavefunction 
vanish when a trial wavefunction does. This gives the 
best upper bound consistent with the assumed nodes. 
With the nodal assumption, the fixed-node approach 
does not have the fermion sign. Since there are sev- 
eral recent reviews on ground state methods (Schmidt 
and Kalos,1987; Senatore and March, 1994; Hammond, 
Lester, and Reynolds, 1994; Anderson, 1995), I will not 
discuss those here. 

Path Integrals have some significant advantages over 
zero temperature methods (as well as disadvantages, of 
course). Among the advantages is the absence of a trial 
wavefunction which means that quantum expectation 
values, including ones not involving the energy, can be 
directly computed. For the expert, the lack of an impor- 
tance function may seem a disadvantage; without it one 
cannot push the simulation in a preferred direction. How- 
ever as the quantum system becomes more complex, it 
becomes increasingly difficult to devise satisfactory trial 
functions. It is my opinion that if quantum simulations 
are to attack the type of complex physical situations in 
which classical simulations routinely deal with, it is bet- 
ter to have a formulation without a trial function. Only 
the Hamiltonian should enter. Of course, the explicit 
formulation at finite temperature also makes comparison 
with experiment more direct. 

I will briefly mention various ways around the fermion 
sign problem, but concentrate on the restricted fermion 
path integral method. This approach has not been de- 
scribed in much detail previously, so I will take the op- 
portunity to do so here. In addition to giving some sim- 
ple examples of the application of restricted paths, I will 
discuss a few practical aspects of fermion path integrals 
such as sampling and nodal actions. Finally, I will briefly 
review some of the applications in order to point out 
the advantages of RPIMC, the future directions, and its 
problems. I shall only discuss continuum fermion systems 
in equilibrium and not discuss lattice models or dynam- 
ical aspects. Gubernatis and deRaedt cover those topics 
elsewhere in this volume. 



However, the straightforward application of those tech- 
niques to Fermi systems means that odd permutations 
subtract from the integrand. This is the "fermion sign 
problem" which I will derive in some detail later on. Path 
integral methods as rigorous and successful as those for 
boson systems are not yet known for fermion systems in 
spite of the activities of many scientists throughout the 
last four decades. On the other hand there has been sub- 
stantial development of zero temperature fermion simu- 
lations using the fixed-node method(Anderson,1976). A 
well-known example is the calculation of the correlation 



II. IMAGINARY TIME PATH INTEGRALS 

In this lecture, we shall consider a system of N nonrel- 
ativistic particles interacting with a potential V(R). The 
Hamiltonian is: 



-H = -J2 X ^'j + V(R) 



(2.1) 



i = l 



where we shall use the definition A,- = S 2 /2m,-. 
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Thermodynamic properties are averages over the ther- 
mal iV-body density matrix which is defined as a thermal 
occupation of the exact eigenstates <f>i(R): 

p(R,B!;fi)=^2rdR)e- pEi M^)- (2-2) 

i 

The partition function is the trace of the density matrix. 

Z(I3) = e-P p = f dRp(R, R;l3) = J2 ( 2 -3) 
J i 

Other thermodynamic averages are obtained as: 

(O) = Z(p)- 1 j dRdR'{R\0\R')p(R',R; 0). (2.4) 

As mentioned above, the density matrix can be calculated 
using path integrals. 

Now let us consider how particle statistics are ex- 
pressed in path integrals. For systems of identical par- 
ticles, the states can be classified into symmetric and 
antisymmetric states. The fermion density matrix is de- 
fined by restricting the sum to be only over antisymmet- 
ric states. (Similarly for other symmetries such as mo- 
mentum or spin.) We shall denote the statistics of the 
particles by subscripts: pp will denote the fermion den- 
sity matrix, ps the boson density matrix, pp the boltz- 
mannon (distinguishable particle) density matrix, and p 
any of the above density matrices. 

Note that for any density matrix the diagonal part is 
always positive: 

p(R,R;0)>O (2.5) 

so that Z~ x p(R, R; (3) is a proper probability distribu- 
tion. It is the diagonal part which we need for many ob- 
servables, so that probabilistic ways of calculating those 
observables are, in principle, possible. 

For the moment, let us consider a single component 
system of spinless fermions. (By spinless fermions, we 
mean that the spatial wavefunction is antisymmetric with 
respect to interchange of all pairs of spatial variables.) 
Let V be one of the N\ permutations of particle labels. 
Then each of the fermion eigenstates has the following 
property: 

4>{VR) = (-l) V HR). (2.6) 

Using Eqs. (2.2) and (2.6) the density matrix has the 
following symmetries: 

p{R,R'-p) = p{R\R-p) 
PF (R, R'- 13) = (-lf PF (VR, R'; p) (2.7) 
Pf(R,R';P) = {-lfp F {R,VR';p). 

One can use the permutation (or relabeling) operator 
to construct the path integral expression for the boson or 



fermion density matrix in terms of the Boltzmann density 
matrix: 

Pb/f(R,R';P) = -^J2(±lfp D (VR,R';l3). (2.8) 
' V 

More generally, one uses some projection operator to 
select a desired set of states from the distinguishable 
particle density matrix which contains all states. From 
Eq. (2.7) we could (anti)symmetrize with respect to the 
first argument, the last argument or both. This connec- 
tion between the boltzmannon density matrix and the 
Bose/Fermion density matrix is important because it is 
the boltzmannon density matrix that is built naturally 
from paths. In this lecture, except for how paths close, 
particles are generally considered to be distinguishable. 
This is in contrast to the second-quantized philosophy, 
where one always works with an antisymmetric basis. 
The distinguishable basis has advantages for numerical 
computation. 

An alternative definition of the density matrix is by its 
evolution in imaginary time, the Bloch equation: 

which obeys the boundary condition at t — for boltz- 
mannon statistics: 

p D {R,R';Q) =8{R-R') (2.10) 

or for Bose or Fermi statistics: 

Pb,f(R, R'; 0) = ^(±1) V S(VR - R'). (2.11) 
' V 

The high temperature boundary condition is an 
(anti)symmetrized delta function. 

Path integrals are constructed using the product prop- 
erty of density matrices: 

p{R ,R'2,Pi+fo) = J dR lP {R a ,R 1 -(3 1 )p{R u R 2 -(3 2 ). 

(2.12) 

The product property holds for any sort of density ma- 
trix. 

If the product property is used M times we can relate 
the density matrix at a temperature /3 _1 to the density 
matrix at M/3" 1 : 

Pb/f(Ro,Rm,P) = 4yl](±l)^ IdRx-.-dRM-i (2.13) 
• <p 

Pd (VRo,Ri;t) . . -Pd(Rm-i,Rm;t). 

The sequence of intermediate points {Ri, R*2, ■ ■ ■ , Rm-i} 
is the path, and the time step is r = /3/M . 
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As the time step gets sufficiently small, we can write 
down an explicit expression for the density matrix pp 
and thereby an explicit expression for p(Ro, Rm', ft), but 
one with lots of intermediate integrals and the permu- 
tational sum to perform. The Trotter theorem tells 
us that for sufficiently small r we can assume that 
the kinetic and potential operators commute so that: 
e _T ^ = e TT e~ TV . Define the (boltzmannon) action as 
S d (R,R';t) = -\ii\p d (R,R';t)]. Then for small time 
step the action is: 



S D (R, R'; r) = ^ ln(4^Ar) + ^— - 



A2 



R') 



1 



-(V(R) + V(R'j). 
(2.14) 



This is known as the primitive approximation to the ac- 
tion. The form of the action is analogous to the poten- 
tial energy of classical "polymer" system with harmonic 
springs between neighboring beads and an interpolymer 
potential between different chains. 

The boson action is real, but it is expressed as a sum 
over permutations. For large N it is not possible to eval- 
uate direcly the sum since it has TV! terms. It is bet- 
ter to leave the bosonic symmetrization as an explicit 
boundary condition on the paths and to sample the per- 
mutations as well as the paths. We will follow the same 
philosophy with restricted fermion paths, the reasons not 
being the difficulty of evaluating the resulting determi- 
nant (that is easier for fermions than for bosons) but 
to avoid the minus signs. For bosons or fermions we can 
(anti)symmetrize anywhere along the path as many times 
as we like. However relabeling is only necessary once. By 
convention we will just relabel the first step. 

The Feynman-Kacs formula for the density matrix re- 
sults from taking the limit M — > oo. The action for the 
path can be interpreted as as Brownian (diffusion) mo- 
tion weighted by potential energy along the path: 



p D {Ro,Rp\P) = p a D{Ra,Rp;P) ( exp 



rP 1 
[- j dtV(R t )} 



BW 



(2.15) 



where the brackets mean to average over all Brownian 
(Gaussian) walks and p° D is the free particle density ma- 
trix (i.e. V - 0). 



III. THE DIRECT FERMION PATH INTEGRAL METHOD 

In the direct fermion method one sums over permuta- 
tions just as for bosonic systems. Odd permutations then 
contribute with a negative weight. The direct method 
has a major problem because of the cancellation of pos- 
itive and negative permutations. This was first noted 
by Feynman and Hibbs (1965), page 292-293 who, after 



describing the path integral theory for boson superfluid 
4 He, noted: 



The (path integral) expression for Fermi 
particles, such as He 3 , is also easily written 
down. However, in the case of liquid He 3 , the 
effect of the potential is very hard to evaluate 
quantitatively in an accurate manner. The 
reason for this is that the contribution of a cy- 
cle to the sum over permutations is either pos- 
itive or negative depending on whether the 
cycle has an odd or even number of atoms in 
its length L. At very low temperature, the 
contributions of cycles such as L = 51 and 
L = 52 are very nearly equal but opposite 
in sign, and therefore they very nearly can- 
cel. It is necessary to compute the difference 
between such terms, and this requires very 
careful calculation of each term separately. It 
is very difficult to sum an alternating series 
of large terms which are decreasing slowly in 
magnitude when a precise analytic formula 
for each term is not available. 

Progress could be made in this problem if 
it were possible to arrange the mathematics 
describing a Fermi system in a way that cor- 
responds to a sum of positive terms. Some 
such schemes have been tried, but the result- 
ing terms appear to be much too hard to eval- 
uate even qualitatively. 



The (explanation) of the superconduct- 
ing state was first answered in a convincing 
way by Bardeen, Cooper, and Schrieffer. The 
path integral approach played no part in their 
analysis, and in fact has never proved useful 
for degenerate Fermi systems, [my italics] 

As Feynman argued above, for the calculation of any 
operator by the direct fermion Path Integral Monte Carlo 
there is a tremendous loss of efficiency over the boson 
case. Here I demonstrate that a distribution having both 
positive and negative regions will have an exponentially 
vanishing (in a sense to be specified) signal/noise ratio 
in a Monte Carlo calculation. 

First, we need to show that the optimal sampling func- 
tion is the absolute value of the integrand which for 
fermion path integrals is just the the bosonic path prob- 
ability distribution. Consider the evaluation of an expec- 
tation value for an operator O: 



(0) = 



JcIxtyO 
J dxir 



(3.1) 



where it is a function with both positive and negative 
pieces and J dx indicates not only an integral over coor- 
dinates, but also the sum over the permutations. Now 
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consider doing importance sampling with some arbitrary 
sampling function P{x). By sampling function, we mean 
that P(x) > and J dxP{x) = 1. The estimator for O 
can be written as: 



Efci w(xi)Q(xi) 



(O) = lim 

M— *oc 



(3.2) 



where the weight is: wixi) = 7r(a;,)/P(a;,) and {xi} is a 
set of M points sampled from P{x). The variance of the 
estimator is: 



vq = 



M[J 



dx 



[*(x)(0(x)-(Q))]< 
P(x) 



(3.3) 



For simplicity, we have assumed the sample points are 
statistically uncorrelated with each other. Otherwise we 
would have to include an additional correlation factor. 

Let us now vary the distribution P{x) to minimize the 
error. The sampling function minimizing the variance is 
found to be proportional to: 



Pq{x) oc \tt(x)(0(x) 



(0))\- 



(3.4) 



The optimal distribution depends on the desired operator 
O. Normally, one wants a sampling function which is 
good for a variety of properties and does not depend on 
the details of a single operator. A physical operator O , 
is usually more smoothly varying than tt(x), which is 
sharply peaked in path space. Also, (O) is not precisely 
known so it would have to be self-consistently determined 
by the sampling. For these reasons, one should ignore 
the second factor and choose P(x) cc |7r(a;)|. With this 
choice the weights are w(x) = sign(ir(x)) = ±1. Then 
substituting into the variance equation we obtain: 



J G?a;|7r(a 



M[fdxir(x)] 2 



j dx\-K{x)\{0{x) - (0)) 2 . (3.5) 



Now for a bosonic system: 7r(a;) > so we will write the 
bosonic variance as: 

v o = ~mJ\A J Hx)l{ ° {x) ~ {0)? - (3 - 6) 

Then we can write the fermion variance as: 

v Q = v B H (3.7) 
where the efficiency (related to the sign of 7r) is 
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"M+ - M_" 


1/hJ 




M 



€= ^ = ""^f^ • (3-8) 



M + j M is the average time the simulation spends in the 
positive region of 7r and M_ / M the average time in the 
negative region. The fermion efficiency is proportional 
to the square of the average sign: the excess of positive 
sampled points over negative sampled points. 



Now let us specialize to the direct fermion method in 
Eq. (2.14). The sign of -k is equal to the sign of the 
permutation, so the efficiency is simply the number of 
even permutations minus number of odd permutations. 
In fact, we can say the efficiency is exactly equal to: 



~Zb~ 



exp[-2 / 3( J F F - F B )] 



(3. 



where Zp and Zb refer to partition function and Fp and 
Fb to the total free energies for Fermi and Bose statistics 
respectively. Of course the free energies are proportional 
to the number of particles. At high temperatures, it can 
be shown (Landau and Lifshitz, 1970, sec. 55) that: 



£ = exp[-2pA r (27rA/3) 3/2 ] 



(3.10) 



At the degeneracy temperature: £ ss exp(— iV)! Loh et 
al (1990) and Newman and Kuki(1992) show figures of 
the average sign for lattice models and for two fermion 
systems. The direct fermion method, while exact, be- 
comes exceedingly inefficient as (5 and TV increase, pre- 
cisely when the physics becomes interesting. 

Sampling the absolute value of an integrand is the best 
we can do without using particular properties of the op- 
erator O. Correlation between the points of the random 
walk used to sample 7r will only further reduce the effi- 
ciency. (In the fermion case one might do a bit better by 
enhancing the probability for tunneling between the pos- 
itive and negative regions of the integrand.) The above 
estimate is the best case of how the variance depends 
on computer time, temperature, and the number of par- 
ticles. One can see a loss of efficiency just due to the 
change of sign of n. Thus, with direct sampling, simula- 
tions of many-fermion systems at low temperatures are 
essentially hopeless because the computational efficiency 
goes to zero too rapidly. 



IV. THE RESTRICTED PATH INTEGRAL METHOD 

We now derive the restricted path identity: that the 
nodes of the exact density matrix determine the rule by 
which one can take only paths with the same sign. For 
the diagonal density matrix we can arrange things so that 
we only get positive contributions. We will do the proof 
for any symmetry, not just antisymmetry, although this 
lecture is only concerned with fermions. 

Suppose pp is the density matrix corresponding to 
some set of quantum numbers which is obtained by using 
the projection operator A on the distinguishable parti- 
cle density matrix. Then it is a solution to the Bloch 
equation (2.9) with the initial condition: 



p F (R, R*;0) = AS(R-R^). 



(4.1) 
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The point R* , the right leg of the density matrix, known 
here as the reference point, will be held fixed. Since A is 
a projection operator we have: pp(R t , R t ; ft) > 0. 

The argument that one can restrict that path integral 
is relatively simple, but for convenience we break it into 
several steps. 

1. The solution to the Bloch equation (2.9) is uniquely 
determined by its boundary conditions. Consider 
some region T of "space-time" . Then the solution 
in the interior is uniquely determined by the values 
on the boundaries of T. More than this, to deter- 
mine the value of the density matrix at a space-time 
point (R, ft), all that are needed are the boundary 
values for ft' < ft. 

A thought experiment with a computer big enough 
to store a fine grid in 3iV-dimensional space gives 
a intuitive proof of this result. Suppose at one 
time slice we know all values of Pf{R) inside the 
spatial domain. Then we can use the Bloch equa- 
tion and the boundary conditions to advance to the 
next time step. Hence the present plus the bound- 
ary conditions determine the future. See Ceperley 
(1991) for the proof of the uniqueness property. 

2. The nodes of Pf{R, R*',t) carve up space-time into 
a finite number of nodal cells, Tk(R t ). Let us define 
a node- avoiding path, as a continuous path R t for 
< t < ft for which p F (R t ,R t ;t) ^ for all < 
t < ft. We say that two points are in the same 
nodal cell if they are connected by a node-avoiding 
path. The collection of all space time points (R, t) 
connected by some node-avoiding path make up a 
nodal cell Tk(R t ). The surfaces of Tk(R t ) consist 
of points such that pp(R, R*;t) = 0. 

From i), we can solve the Bloch equation inside 
each nodal cell separately by specifying the initial 
conditions at t = from Eq. (4.1) and zero bound- 
ary conditions on the surface of the nodal cell. We 
refer to the union of all nodal cells as T(ii«). 

3. Now let us reconstruct the path integral expression 
for the density matrix taking into account the zero 
boundary conditions. To enforce the zero boundary 
condition, we can insert an infinite repulsive poten- 
tial precisely at the nodal surfaces. But that will 
eliminate the contribution of any walks which hit 
or cross the node. Hence, demanding that a density 
matrix vanish on a surface is equivalent to restrict- 
ing the path integral to be only over node-avoiding 
walks. 

Thus we have proved the Restricted Path Integral iden- 
tity: 

p F {Rp,R t ;ft) = I dR o p F (R o ,R*;0) 



i dR t e~ s[R ^ (4.2) 

JRo->RpeT:(R,) 

where the subscript means that we restrict the path in- 
tegration to paths starting at Ro, ending at Rp and are 
node-avoiding (those for which p F (R t , R*;t) ^ for all 
< t < ft.) The weight of the walk is p F (R , R*; 0). It 
is clear that the contribution of all the paths for a single 
element of the density matrix will be of the same sign; 
positive if pp(Ro, R*;0) > 0, negative otherwise. This is 
what we set out to show. In particular, on the diagonal 
all contributions must be positive. 

Important in this argument is that the random walk 
is a continuous process (the trajectory is continuous) so 
we can say definitively that if sign of the density matrix 
changed, it had to have crossed the node at some point. 
This presents a technical problem for discrete time paths 
since we must decide whether the path crosses and crosses 
back in between the sampled times. We take up this 
issue later. Lattice models or non-local Hamiltonians do 
not have continuous trajectories so this method is not as 
straightforward for those models. 

The restricted path identity is one solution to Feyn- 
man's task of rearranging terms to keep only positive 
contributing paths for diagonal expectation values. The 
"bosonic" path integral formulation can be applied to 
fermion path integrals; all that is required is to take a 
subset of the bosonic paths. In principle, there exists 
a way to solve the "sign problem" ! We shall see that 
it is important to allow non-trivial, even permutations. 
Macroscopic even permutations are directly related to 
Fermi liquid behavior. 

The problem we now face is that the unknown density 
matrix appears both on the left-hand side and on the 
right-hand side of Eq. (4.2) since it is used to define the 
criterion of node-avoiding paths. To apply the formula 
directly, we would somehow have to self-consistently de- 
termine the density matrix. In practice what we need to 
do is make an ansatz, which we call pj* , for the nodes of 
the density matrix needed for the restriction. The trial 
density matrix is used to define trial nodal cells: (R*)- 
Using a trial nodes we generate a better approximation 
to the density matrix using Eq. (4.2) with the trial re- 
striction: 

Pr(Rp,R.;p)= j dR p F (R , R,; 0) 

I dR t e- s ^ Rt \ (4.3) 

Hence, p^ (R! , R; ft) is a solution to the Bloch equation 
inside the trial nodal cells, and it obeys the correct ini- 
tial conditions. It is not an exact solution to the Bloch 
equation (unless the nodes of p^ are correct) because it 
has possible gradient discontinuities at the trial nodal 
surfaces . 
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By multiplying together p? and pp, integrating over 
time and a common spatial variable, and using their re- 
spective definitions, one can show that p? is related to 
the exact solution by the following integral equation: 

p T (R u R 2 ;t)=p F (R u R 2 ;t) 

+ dt' J dR'p F {Ri,R'-t-t')K{R' ',R 2 ;t') 

where the kernel is non-zero only on the nodes of p? 
where it equals the discontinuous derivative of p? across 
the node. 

K(R 1 ,R 2 ,t) = XS(p T (R 1 ,R 2 ;t)) 

V\p T {Rt,R 2 ;t)-p T {R^,R 2 ,t)] (4.5) 

Before we discuss forms for the trial density matrix, 
let us return to the more general case where one has spin 
h fermions. Assume that there is no magnetic field so 
that the Hamiltonian is independent of spin. Then we 
can quantize the spin axis in the z direction. Suppose 
we want to calculate the partition function Z m in an 
ensemble where S z is fixed to be m. It is not hard to see 
that this equals: 

S N {a,Va)S{Y J ^i,m) J dRp D {R,VR)^ (4.6) 

Now simply relabel the coordinate indices in the integra- 
tion so that particles 1 through iVf = N/2 + m have up 
spin and particles + 1 through N have down spin. 
Then clearly all arrangements of the spin variable in the 
trace contribution the same amount to the partition func- 
tion and we can just calculate one term and multiply 
by the number of ways of assigning spins. There are 
7V!/(iV-[-!iV|!) such arrangements. Define the following 
operator to antisymmetrize over the up spins and the 
down spins individually. 

where V j' (4-) operates on the up (down) spin coordi- 
nates only. With this one can calculate spin-independent 
correlation functions in the fixed S z ensemble. 

If we want to calculate the total partition function: 
Z = Z m the magnetization must be sampled. To 
average over all spins one could occasionally attempt a 
spin flip from one fixed m ensemble to ra±l and accept 
or reject that flip, depending on whether the path is now 
legal. 

We can also do a fixed S ensemble by using a differ- 
ent projection operator. For example one can use the 
relationship: 



Z s = (2S+l)(Z m -Z m+1 ) (4.8) 

and try to calculate the differences of partition functions. 
Equivalently one could use a projection operator for the 
total spin. This is discussed in detail in the recent work 
of Lyubartsev and Vorontsov-Velyaminov (1993). As an 
example, suppose we want to calculate the partition func- 
restricted to the spin-singlet states of a four electron 
system like the Be atom. Then one wants to sum only 
over 12 of the total number of 4!=24 permutations. The 
identity permutation gets a relative weight of 2, the three 
possible double pair exchanges such as (2,1,4,3) get a rel- 
ative weight of 2, and the eight possible triplet exchanges 
such as (2,3,1,4) get a relative weight of -1. These weights 
specify the projection operator that is used for the ini- 
tial value of the density matrix. What is not clear to 
me is how these approaches will scale in the number of 
fermions. 

V. THE REFERENCE POINT 

We call R* the reference point and it plays a very spe- 
cial role in restricted path integrals since it is the value 
of the density matrix with respect to the reference point 
that restricts the paths. Averages such as the density 
can only be taken at the reference point. This is different 
from boson or distinguishable path integrals; for them all 
time slices are equal by an obvious time-slice symmetry, 
and all time slices are equivalent in taking averages. If 
we write expectation values in terms of derivatives of the 
partition function (e.g. the kinetic energy is the mass 
derivative of the partition function), one recovers an ex- 
pectation value involving all the time slices, not just one 
containing the reference point. However, evaluating the 
derivatives can be quite complicated. 

By a "time-independent" or "ground-state" restriction 
is meant that the restriction does not depend on the ref- 
erence point. This is achieved by using an antisymmetric 
trial wavefunction ^^(R) and requiring that if>T(R t ) ^ 
throughout the path. This is identical to the ground state 
fixed-node method (diffusion or Green's Function Monte 
Carlo). The algorithm is considerably simpler than us- 
ing time-dependent nodes and time-slice symmetry is re- 
stored. However, calculations with ground-state nodes on 
liquid 3 He gave a poorer description of the properties of 
the liquid at non-zero temperature (Ceperley, 1992). It is 
not clear in what sense this time-independent restriction 
would give proper finite-temperature effects, even though 
one is using path integrals. I speculate that any exact 
fermion description at finite temperature must break the 
time-slice symmetry if it is to solve the sign problem. 

There is a simple trick for using two reference points 
instead of one. This improves the method in several 
ways. Let us square Eq. (4.2). This means we con- 
struct two paths, one from Rq to Rp and another from 
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R' to Rp both in time f3. The weight of the product of 
the two paths will be: pp(Ro, R f ;0)pp(R' o , R f ;0). When 
we have integrated over all paths, and the two initial 
positions Rq,R' q and the reference point, we will have 
determined the density matrix J <IR*Pt(R*, Rp\ P) 2 ^ 
PF{Rp, Rp',^P), that is at a temperature half as big. 
Note that both of the paths are node avoiding with re- 
spect to the same reference point, R* and end up at 
the same point Rp at the same time. Hence, the sign 
of pf{Rq, R*]0) equals the sign of pp(R f , Rp, /?), which 
equals the sign of pp(R' Q , R*;0). The contribution of 
the double path is positive. Also note a point VR* 
has the same effect on the restriction as does R* since 
Pf(R, R*) = implies that pf(R, VR*) = 0. This means 
we can work with only a single projection and use a 
weight for the double path since: 

p F (Ro,R' o ;0) = J dR*pF(R o ,R*;0)p F (R' o ,R*;0). (5.1) 

The net effect of doubling the time is to redefine how 
we measure time differences in the restriction. Let the 
reference point be defined as having "time" value of zero. 
Then the time argument that we should use for the re- 
striction (the restricted time) is: 

(t for < t < (3/2 
r ~\(3-t for /3/2<t</3 [ ° ' 

(Now we are using an inverse temperature (3 not 2/3.) 

Time doubling is an improvement because if we have 
accurate nodes down to a temperature T, we can do 
accurate simulations down to T/2. Also, we can use 
both Ro and Rp/2 to calculate expectation values. Fi- 
nally, we have restored (imaginary) "time-reversal" in- 
variance. If our nodes are incorrect, then in general: 
Pt(Ri,Ro',P) 7^ Pt {Ro, Ri, P)l the symmetry of Eq. 
(2.7) could be violated. By using two paths from Rq 
to -R/3/2 we are performing the integral: 

p T {R 1 ,R 2 -(3) = J dR*p T {R u R*-pl2)p T {R 2 ,R*-pl2) 

(5.3) 

which is clearly symmetric in R\ and R%. 

The time doubling cannot be repeated without reintro- 
ducing the sign problem. Also, one gets a sign problem 
if the two time arguments that are multiplied together 
are different, so this method cannot be used to determine 
imaginary time correlation functions without having sign 
problems. But it is well worth doing it once. 

VI. AN EXAMPLE OF RESTRICTED PATHS 

To get a feeling for restricted paths let us consider 
the problem of molecular hydrogen. We will work in 



the Born-Oppenheimer approximation so the two pro- 
tons in a single hydrogen molecule are represented by 
two spin 1/2 particles interacting with an attractive po- 
tential. The total spin (S) and total orbital angular mo- 
mentum (L) are good quantum numbers. The spin 
state is called para-hydrogen, and must have an even 
value of L to keep the molecular wavefunction antisym- 
metric in spin and coordinates. The spin 1 states are 
called ortho-hydrogen and they must have odd values of 
L. Often the hydrogen cannot easily change its spin state, 
so that para- and ortho-hydrogen can be considered as 
separate chemical species, for a time at least. They can 
change their angular momentum values with collisions 
with other molecules but not easily their spin. A third 
possibility of statistics is if the two nuclei are different 
particles, e. g. a proton and a deuteron, in which case 
they obey distinguishable particle or Boltzmann statis- 
tics. 





a 




• 






b 










c 




• 





FIG. 2. Three types of paths that contribute to the partition 
function of a single hydrogen molecule. Shown is the vector 
separation between the two atoms; the sphere represents the 
repulsion when the two atoms approach each other. The two 
small circles, represent the reference point, the origin of time. 
The irregular curve represents the path; it starts and ends at 
a reference point. The horizontal line represents the nodal 
plane in the ortho-hydrogen state. 
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Let us suppose that the particles are massive enough 
that the relative coordinate is almost fixed at a given 
radius ro « 0.75 A. Hence the relative coordinate r = 
ri — r2 is almost fixed on the surface of a sphere. In 
figure (2) we show four types of paths. The reference 
point is the large dot. 

Now consider what we need to do to calculate the parti- 
tion function for the three types of statistics. Distinguish- 
able particles are the simplest: allow all paths returning 
to the starting point (type A in the figure). For ortho- 
and para-hydrogen, we can use parity to project out the 
correct states. This generates paths of type B which 
end up at the opposite pole in relative coordinates. For 
para-hydrogen the direct method would be to sum over 
all paths of types A and B. Ortho-hydrogen would be to 
sum over paths of type A but subtract the contribution 
from paths of type B. 

In the case of ortho-hydrogen it is easy to calculate 
the exact nodes of the density matrix. In relative coor- 
dinates any wavefunction has the angular factor Y; m (r) . 
Then in the sum over the quantum states for different 
m in Eq. (2.2), using the addition formula for spherical 
harmonics we obtain a factor P;(r -r t ). Since all the odd 
Legendre polynomials vanish when their arguments van- 
ish the ortho-density matrix vanishes when r ■ r* — 0. In 
general, additional nodes would be possible, but the hy- 
drogen molecule has only this one planar node. Paths of 
type C are node-crossing as opposed to the node-avoiding 
path A. 

To summarize, we add the following classes of path for 
the different statistics of the hydrogen molecule: 

1. Distinguishable hydrogen: A+C 

2. Para-hydrogen: A+B+C 

3. Ortho-hydrogen (direct method): A-B+C 

4. Ortho-hydrogen (restricted method): A only (r t • 
r, > 0) 

The reason that restricted path integrals give the same 
value is that paths of type B and C can be paired together 
and canceled off against each other. This is because the 
flux of paths is the gradient of the density matrix at 
the node and since the gradient is continuous across the 
node, the positive paths crossing at a given nodal point 
will precisely cancel against the negative paths. 

Hence the restricted paths are limited to be in a half- 
space. Note that there is no definite location of this half- 
space. Its position depends on the reference point, even 
at zero temperature. This is because the ground state 
of S=l is three-fold degenerate. Isotropy is restored by 
averaging over the reference point position. 

Consider now a crystal made of hydrogen molecules. 
To a good approximation one can neglect the angu- 
lar coupling of the spins in constructing the nodal sur- 



face. From the point of view of restricted path in- 
tegrals, ortho-hydrogen is more easily orientable than 
para-hydrogen because its restricted paths are naturally 
dumb-bell shaped. Runge et al (1992) have used this 
approach in PIMC simulations of both ortho- and para- 
hydrogen. On the other hand, Buch (1994) in his PIMC 
calculations, explicitly diagonalized the rotational Hamil- 
tonian at each step, a considerably more complicated and 
more approximate way of handling ortho-hydrogen. 

VII. NODES OF THE DENSITY MATRIX 

The only uncontrolled approximation in the Restricted 
Path Integral Method is the restriction, the rule by which 
we allow paths. Clearly the success of the method hinges 
on the choice of this restriction. The situation is not very 
different from that of classical Monte Carlo or Molecular 
Dynamics simulations. Except for a few isolated cases, 
the true classical Hamiltonian describing a real physical 
system is not very well known (say compared to a typi- 
cal thermal energy of 100 K ss 0.01 eV ss 0.35 mH). I am 
not referring to simplified models like the Ising model 
or a Lennard- Jones liquid but about water or a poly- 
mer. In principle, the classical dynamics is controlled 
by the true Born-Oppenheimer potential. In practice, a 
very simplified potential is used, either a semi-empirical 
one or one from on-the-fly LDA calculation as Car has 
described. One knows that the true trajectory deviates 
very quickly from the calculated trajectory and that de- 
tailed properties such as the melting point depend on the 
actual potential. However, simulations are still useful be- 
cause they are true simulations of a nearby model system. 
Hopefully, the model is in the right universality class. 
Simulations of good models are useful because they are 
the only way to understand reliably many-body effects. 
Even if some of the quantitative details are inaccurate, if 
we can characterize the nodal restriction sufficiently well, 
the simulations will be useful in understanding nearby 
strongly-interacting fermion systems. 

Let me explore further this analogy between classi- 
cal simulations and RPIMC. The dimensionality of the 
unknown "potentials" is not so different. The Born- 
Oppenheimer classical potential is a real function of 3iV 
variables. The nodal restriction is a Boolean function (i.e. 
only the sign matters) of 6N + 1 variables (counting the 
two legs of the density matrix and imaginary time). An- 
other similarity is that at finite temperature, the nodal 
surface is local. By this we mean that the nodal surface 
is determined by the equality of even and odd local per- 
mutations; local in that each particle can only move on 
the order of a thermal wavelength. For classical neutral 
atoms with localized electrons, one generally assumes the 
potential is a pairwise sum with smaller and higher three- 
body contributions, so the many-body potential can be 
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built up from the properties of small clusters. In both 
cases, permutational and other symmetries can be used 
to reduce drastically the phase space where we need to 
determine the "potentials." 

The nodal surfaces have hardly been investigated, so 
there is much to be done here. It is my opinion that 
we should not be put off by the RPIMC method just 
because we do not know the exact restriction, but we 
should take inspiration from what has happened in clas- 
sical simulations and push ahead with a nodal ansatz, 
all the while trying to improve our knowledge of fermion 
nodes. Hopefully, the trial nodes are good enough to 
put the model in the right universality class. What we 
don't have is the century of accumulated wisdom about 
what nodal properties are important. In the rest of this 
section, we review what is known about the nodes, and 
various approximations that could be used. 

In an actual calculation, one does not make a geometric 
interpretation of the nodes as we did in the two-particle 
example. Instead one computes a trial density matrix 
and uses that to decide whether a given path is to be 
allowed. Paths for which the trial density matrix are 
negative or too close to zero are rejected. I will discuss 
this aspect in the next section. How can we choose a fully 
antisymmetric trial density matrix which can be quickly 
computed during the Monte Carlo random walk? Using 
the projection, the fermion density matrix is: 

PF (R, r, ■ t ) = Ap D (R, R* ; t ) c< (- 1 ) V pd (PR, R* ; t ) . 

V 

(7.1) 

Only the sign matters so we can drop any constants. This 
formula is exact but not useful unless a) we have a good 
expression for pp down to temperatures as low as 2T and 
b) we can evaluate quickly the permutational sum. Effec- 
tively, the only convenient antisymmetric density matri- 
ces for the many-body case are determinants of one form 
or the other. 

For the moment, let us consider the reference point R* 
and the inverse temperature t as fixed parameters. Then 
the nodes have dimension 3N — 1 since a single equa- 
tion specifies whether R t is on a node. One property 
that holds true in general, is that when two fermions 
have both the same spin and the same spatial coordi- 
nates, all the wave functions and hence the density ma- 
trix must vanish. Hence, for any pair of fermions with 
the same spin, the hyperplane defined by the equation: 
r, = Tj is on the node. Since these are three equations 
(in three dimensions) the "coincident hyper-planes" have 
dimensionality 3N — 3. The coincident planes are fixed 
hyper-points lying on the nodal surfaces which have a di- 
mensionality two larger. For quantum mechanics in one 
dimension, the coincident points exhaust the nodal sur- 
faces so that one knows the exact restriction. This has 



misled many researchers into believing that the situation 
in higher dimension was not much more difficult. Such 
is not the case, at least as far as we know. For fermions 
in two or three dimensions, symmetry is not sufficient 
to determine the position of the nodes. Their position 
depends in a non-trivial way on the potential. (There is 
also a significant difference in complexity between 1 and 
higher dimensions. In ID the sign of the density matrix 
is equivalent to ordering on a line, so computing the re- 
striction is a completely local evaluation, while in higher 
dimension computing even the free particle density ma- 
trix costs C[iV 3 ] operations.) Knowing the coincident 
points are on the nodes means that if one particle is dis- 
placed roughly an interparticle spacing, typically one will 
cross a node. Later, we will see the connection between 
the nodes and the momentum distribution. 

Now, let us discuss the nodal surfaces of non- 
interacting fermions. Let v e (r) be a single-particle exter- 
nal potential. The distinguishable particle density matrix 
is then a product of solutions of the single-particle Bloch 
equation: 



dg(r,t) 
dt 



[-AV 2 + . e (r)]j(r,t) 



with the boundary condition: 

#(r,r»;0) = £(r-r»). 



(7.2) 



(7.3) 



Then using the antisymmetric projection operator and 
the definition of a determinant, we find for the spinless 
case: 

Pf(R, R*',t) = ^ydettfffc.rj,,;*]. (7.4) 

In the case where the external potential is zero (or a con- 
stant), the single-particle density matrix is a Gaussian. 



#(r, r»;t) oc exp 



(r-r,) 2 



4Xt 



(7.5) 



Note that in periodic boundary conditions, the rela- 
tive distances should be computed using the minimum 
image convention. But if the thermal wavelength y/2Xt 
approaches the size of the box, which can occur at low 
temperatures, additional images may be required. The 
exact single-particle density matrix is a theta function in 
periodic boundary conditions. Only if the exact density 
matrix is used in the determinant, will the many-body 
density matrix approach the Slater determinant of plane 
waves appropriate to the ground state of a uniform liquid. 
(This is an issue for restricted fermion path integrals but 
not for boson path integrals because the "time" argument 
ranges up to 1/(2T) instead of just up to r = lj(MT).) 

At high temperature these nodes have a geometrical 
construction. Consider sitting at a point R and taking 
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the limit as t — > + . For a typical point, a single permu- 
tation V R* will be closest to R and all the other permu- 
tations can be neglected. If V is even, pp will be positive, 
otherwise it will be negative. A point R is on a node if 
it is precisely equidistant from an even permutation and 
an odd permutation. 

(R - V e R*f = (R - V R*f (7.6) 

or 

R ■ (V e - V )R* = 0. (7.7) 

Hence at high temperatures the nodal surfaces are hyper- 
planes. A nodal region is the set of points closer to R* 
than to any other odd permutations VR*. In this sense, 
the nodes are as far away from the reference point as an 
antisymmetric function will allow. 




CONTOUR FROM TO BY OxlO ! 



FIG. 3. A 2D cross section of the density matrix of 161 free 
spinless fermions in a periodic square. The open circles are 
the fixed positions of 160 fermions;the last fermion is free to 
roam around the squares. The contour lines are the nodes 
of the density matrix. The curves respresent the nodes for 
temperatures in the range 2.5 Ef < T < \0Ef where Ef is 
the Fermi energy. 

On the other hand, at low temperatures the free- 
particle nodes are determined by minimizing the kinetic 
energy. They must be as smooth as possible while keep- 
ing the coincident points on them. Hence the nodes 



change character as the evolve from high temperature 
to low temperature as seen in Fig. (3-4). 

I have investigated (Ceperley, 1991) the number of 
different nodal regions into which the full configuration 
space is divided. I have found that the free-particle 
spinless-fermion density matrix in periodic boundary 
conditions divides the total phase space into only a sin- 
gle positive region and a single negative region (except in 
one spatial dimension where the nodes divide the phase 
space into N\ regions.) This means that the restricted 
path partition function includes contributions from all 
N\/2 even permutations. 




CONTOUR FROM TO BY 100 



FIG. 4. A 2D cross section of the density matrix as in the 
previous figure. In this figure the temperature is close to the 
ground state; the temperature ranges from 0.04 Ef < T < 
0.25E F . 

To go beyond free particle nodes, we can first work 
within a one-particle self-consistent approach, such as 
Hartree-Fock. This has the virtues of a) giving a Slater 
determinant of single particle density matrices which will 
be easy to evaluate and b) giving good nodal surfaces at 
zero temperature. It has been found that nodes com- 
ing from Hartree-Fock calculations and LDA calcula- 
tions work quite well for electronic systems. All that 
is required in the above equation (7.2) is to use an ef- 
fective external potential instead of the bare potential. 
It is possible to have the effective potential be refer- 
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ence point dependent i>(r, r»;£) and for the mass be 
an effective mass which is "time-dependent" A(i). This 
roughly corresponds to unrestricted Hartree-Fock, or to 
self-interaction corrected density functional theory. 

Another approach is to take what we know about good 
distinguishable particle density matrices and directly ap- 
ply formula Eq. (7.1). The pair product density matrix 
is: 

Pd{R, R*;t) — p D o{R, R*;t)exp[- ^ u(r ij: r ijtf ; t)] 

(7.8) 

where the pair action u is defined either in terms of the 
exact action of a pair of particles or as the first cumulant 
action. This form is quite accurate for systems interact- 
ing with pair potentials (Ceperley, 1995) at least above 
the degeneracy temperature. 

In the simplest approximation we make the end- 
point approximation for the pair action: u(r, r*;£) = 
h( u d( r ',t) + u d( r *',tj) + C>[i 2 ]- But since the total pair 
action is symmetric under permutations, the pair terms 
will factor out of the antisymmetrization and hence will 
not affect the nodes of the free particle density matrix. 
This explains why the nodes of the free-particle density 
matrix, even though it is over-simplified are a reasonably 
good approximation down to the degeneracy tempera- 
ture. The off-diagonal corrections, which could push the 
nodes away from the free-particle positions, are C[i 2 ], 
while the dominant kinetic energy term is 0[i -1 ]. The 
temperature has to be fairly low for contributions com- 
ing from the potential to matter. This is supported by a 
perturbation calculation of Hall and Price (1991). 

An important feature of accurate actions is that they 
contain off-diagonal pair terms, i.e. the pair action de- 
pends on both R and R* in a more complicated way than 
this. But without the end-point approximation, the re- 
sulting antisymmetric function is not a determinant, so 
that can only be computed directly for systems with a few 
fermions where one can explicitly sum up the permuta- 
tions. One pathological feature of free-particle nodes is 
that there is no coupling in the nodes with unlike spins. 
The nodal surfaces are a Cartesian product of up-spin 
nodes with down-spin nodes. This feature is incorrect 
and possibly very significant in a strongly-coupled sys- 
tem such as liquid 3 He where collective spin excitations 
are important. [A good homework problem is to verify 
that in a four electron system such as the Be atom, the 
restricted paths allow not only the identity permutation, 
the only one allowed in Hartree-Fock, but also the double 
pair permutation (2,1,4,3).] 

In general, we will have to include the effects of the 
off-diagonal terms into the diagonal ones. The principle 
physical effect of the potential is to slow the free-particle 
diffusion. One can define an effective (inverse) mass by: 



y(t) = {(r(t)-r t f)/(6t) (7.9) 

where the average is over distinguishable particle paths 
starting from the point r, . This effective mass is com- 
puted in Ceperley (1995) for liquid 4 He as a function of 
"time" . At short time one has the bare mass but at larger 
times (when the thermal wavelength equals the interpar- 
ticle spacing) the effect of the interparticle interaction 
is to slow down the diffusion by a factor of two. This 
will not affect the nodes at zero temperature, but it will 
rescale the temperature at which those nodes are used. 
In a localized phase such as a solid helium, the effective 
mass diverges so that the single particle density matrix 
remains localized even at zero temperature. 

Let us take a least-squares function-fitting approach 
and say the best single particle function to represent the 
many-body density matrix is the one which minimizes 
the mean squared action: 

X 2 = J dRdR*n(R, R*;t)[S D (R, R*;t) - S A (R, R*;t)] 2 . 

(7.10) 

Here n{R, R*',t) is an measure for the end-points of a 
path chosen to generate a good "sampling" of configu- 
ration space. For example, one can take pairs of points 
{R, R t ) from a RPIMC simulation of the system in ques- 
tions and relabel the R from the distribution of permuta- 
tions pd{V R, R*',i)- In the above equation, Sd(R, R*\i) 
is the exact distinguishable action 

1 N 
S A (R,R,;t) = -[U d {R-,t) + U d (R.-,t)] + Y l *(?i 

i = l 

(7.11) 

is the general form of an approximate action which maps 
into a determinant after antisymmetrization and Ud(R] t) 
is the exact diagonal action. Minimizing the rms action 
with respect to the function s(r, r*;t), we find the opti- 
mal single particle density matrix is proportional to : 

N 

-X)«( r i- r i')>r 1 ,r 1 . (7.12) 

i=2 

where (. . .}r 1 ,r 1 , means an average over [i with (ri,ri») 
held fixed. This equation can be solved by iteration: 
assume the free-particle value for s(r,r*;£) = (r — 
r„) 2 /(4Ai) evaluate the left-hand-side, and iterate. This 
will reduce to the effective mass formula above if we as- 
sume the single particle density matrix is Gaussian and 
we assume the action is normally distributed. 

A more complex way of changing the nodes is to put 
in backflow effects. This has been found to be very suc- 
cessful for liquid 3 He (Panoff and Carlson, 1989) and in 
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the 2D electron gas (Kwon, Martin and Ceperley, 1993) 
getting on the order of 90% of the energy missed by the 
free-particle nodes. To compute the Slater determinant, 
one first transforms the coordinates of all the particles 
on the path to "quasi-particle" coordinates. The sim- 
plest transformation allowed by translation symmetry is: 



Si = r* + ^(r-i - rj^dr; - Tj \;t) 



(7.13) 



where rj(r;i) is the "backflow" function. Then the trial 
density matrix is the usual Slater determinant but writ- 
ten in quasi-coordinates: 



p T (R,R*;t) = det[g(si,Sj*;t)] 



(7.14) 



At high temperature one can show from the Feynman- 
Kacs formula that the backflow function will depend on 
the potential as: 



lim n(r, t) — 
t->o v ' 



\t 2 dv(r) 
2rdr 



(7.15) 



Hence at high-temperature, the classical force pushes the 
nodal surface by an amount proportional to t 2 . At low 
temperature the backflow function rj(r) saturates to its 
zero-temperature form. The difficulty with using a back- 
flow trial density matrix is that all the coordinates get 
coupled together, so that updates of the determinant are 
much more expensive. This is the reason backflow nodes 
have not been often used, even in the ground state. In 
addition, there has been little investigation of ways of 
obtaining the backflow function. 

As we mentioned above in the two particle example, for 
the exact nodal surface the positive and negative paths 
push the nodes equally hard from the two sides. There 
is a sort of equality of nodal surface pressure across the 
nodal membrane. Unfortunately this does not lead to 
a practical general way of constructing nodes because 
the membrane is defined in such high dimensionality. It 
would take too many paths to do that. 

We must bear in mind that all this talk about nodal 
approximation should not make us lose sight of the fact 
that correlation caused by the potential is being exactly 
treated with the distinguishable path integrals. It is only 
the shape of the playing field, not its size or the playing 
rules, that we that we must approximate. Ground state 
calculations of such electronic systems as small molecules, 
the electron gas or many-body hydrogen, show that the 
fixed-node calculation using Hartree-Fock nodes have an 
absolute accuracy of better than 0.01 eV/molecule. Rel- 
ative accuracies are expected to be significantly better. 
One might also suspect that the restricted path errors 
would increase with decreasing temperature, so that the 
ground state would be be worst case. 



VIII. SOME TECHNICAL DETAILS OF RPIMC 
A. The action 

The action for restricted fermions, which we write as 
Sp{Rt, Rt-r', t; t r , R*), is more complicated and depends 
on more variables than the distinguishable action Sd be- 
cause of the lack of time symmetry in the paths. It 
is a function of R* and t in addition to its usual de- 
pendence on R t ,Rt+ T and r. Usually we will not write 
these extra variables, taking them to be understood. For 
ground state nodes (no time dependence) one reverts to 
the bosonic situation; the action does not depend on the 
reference point or reference time. The Feynman-Kacs for- 
mula, Eq. (2.15), for the action arising from the potential 
and the restriction is: 



-U p (R t ,R t+r ) 



exp 



^ dt'V{R t .) 



(8.1) 



The averaging is over all Brownian walks starting at R t 
and ending at Rt+ T and keeping inside T(i? t ). 

The ■primitive-nodal action makes the approximation of 
checking the restriction only at the M sampled points on 
the path. We only check to see whether pr(Rj , R*',tj) < 
for < j < M. To illustrate the problem with the 
primitive nodal action, consider a single particle in a ID 
well of size a: 



V(x) 



-{ 



for < x < a 
+00 otherwise 



(8.2) 



In this model, the nodes are fixed, as they would be for 
any fermion problem at low temperatures. The error of 
the ground state energy, (/? — oo) is shown in figure 5 
versus the time step. In the primitive approximation the 
error decreases as t 1 / 2 . This dependence is easy to un- 
derstand. The energy of a box of size a is A(7r/a) 2 . Using 
the primitive approximation effectively increases the size 
of the box by an amount proportional to the thermal de- 
Broglie wavelength of a time step: (At) 1 / 2 . Shown in 
figure 6 is the density computed with the primitive ac- 
tion. The primitive-nodal action allows the density to be 
non-zero at the nodes. Hence, the energy is decreased by 
a relative amount (At) 1 / 2 /*!. In boson problems, we can 
find actions for which the error in the energy converges 
much faster, as t 3 , so that use of the primitive approxi- 
mation will seriously degrade the overall performance of 
the PIMC algorithm. 

Luckily, one can do considerably better, so that fewer 
time slices are needed to accurately represent the path. 
The action picks up a contribution from the nodes be- 
cause walks can wander back and forth across the nodes 
even though they happen to be on the correct side at the 
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sampled points, R t and Rt+ T - Improved fixed-node sam- 
pling methods have been developed for the ground state 
simulations by Anderson (1976). Here we adapt this to 
path integrals. To systematically improve the action, we 
have to make a sequence of three approximations. 
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FIG. 5. The energy error of a particle in a square well with 
a harmonic potential of size 7r at zero temperature computed 
with four different forms for the action as a function of time 
step. The actions differ in their estimate of the distance d 
of the particle to the wall. The top dotted curve is for the 
primitive action, the distance is set to infinity, the error de- 
creases as t° 5 ; the next lower short dashed curve uses the 
newton estimate with the nodal distance estimated from the 
exact ground state wavefunction so that d = tan(a;), its error 
decreases as r 1 ' 5 . In the next lower dashed curve, the exact 
nodal distance is used but no account is made for the coupling 
between the bosonic action and the nodal action. In the lowest 
dashed curve the distance is estimated as d = min(x,-K — x), 
and there is no bosonic action, the error is due to neglect of 
the mutiple images goes as exp(— ir 2 /4r) 

In first approximation to the average of Eq. (8.1), we 
assume that the potential energy is uncorrelated with 
crossing and re-crossing of the nodal surface. Then the 
action factorizes 

exp[-U F (R t ,R t+T )} = exp[-U(R t ,R t+T ) - U N (R t ,R t+T )} 

(8.3) 

where U is the purely boltzmannon action and t/jv is the 
free particle action due to just the restriction. The sec- 
ond to lowest curve in fig. (5) shows the effect of this 
approximation on the ground state energy for a particle 
in a well with an additional harmonic potential. The er- 



ror in the energy caused by the factorization goes as r 2 . 
The boltzmannon action U has been extensively investi- 
gated for use in the path integrals (Ceperley, 1995) so we 
need say no more about it here. 

The exact nodal action for a particle in a box is easy 
to calculate. In the case of a particle confined only to 
the half-space x > 0, one can solve the Bloch equation 
by the method of images. The total density matrix is 
the difference: pix, x') — p(x, —x') since the difference 
satisfies the Bloch equation and the boundary condition 
at t = and x — 0. Using the form for the free particle 
action the nodal action is: 

U N (x t ,x t+T ) = -Hl-exp(-^T ±L )] (8-4) 

AT 

where dt is the distance of Xt to the node at time t and 
dt+r is the distance of Xt+ T from the node at time t + r. 
Hence the action diverges logarithmically near the node 
and is significantly repulsive in a region on the order of 
\Z\t. See figure 7. 
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FIG. 6. The single particle density in a square well at zero 
temperature with size ir computed with two forms for the 
action for a time step of tE = 0.078 (the corresponding de 
Broglie wave length 0.4). The solid line is the exact density 
sin(a;) 2 . The dotted line is the primitive approximation, its 
density is not zero at the edge of the box. The dashed line 
was computed using the image action. The corresponding 
energies are shown if Fig. 5. 

In particle-in-the-box model, there are boundaries 
(nodes) both at x = and at x = a. Hence, the above 
action is not exactly correct. As in similar problems in 
electrostatics, one must expand in an infinite series of 
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positive and negative images. In the case of the density 
matrix, the images decay as exp(— Z 2 /(4Ar)), so that at 
small r the random walk sees only the nearest wall. If we 
make the approximation of only using the nearest image, 
we get the error shown as the lowest line in figure 5. Its 
error converges exponentially in r. The dominant error 
at small time step is the neglect of the potential, not the 
neglect of the other walls. 

Now, let us return to the many-body case. The free- 
particle nodal action is difficult to evaluate (at least I 
do not know how to do it), so we need to make further 
approximations. We shall assume that locally the nodal 
surface is constant (as a function of t) and is a single 
hyper-plane so we can use the image approximation dis- 
cussed above. 




FIG. 7. The action (solid line) and the energy (rdU/dr), 
dashed line, for a ID particle next to a wall at the origin. 

The remaining non-trivial problem is how to estimate 
the distance d to the node. In the many-body time- 
dependent case we define the distance to the node as 
dt = min(|i? n — R t \) where R n varies over all points 
with Pf{Ru, R*',t) = 0. It is easy to determine an up- 
per bound to dt- If any two particles are moved directly 
towards each other, they will be on a node when they 
coincide. It follows that dt < r,j t /V^. However the 
distance may be much smaller than this. To get an esti- 
mate, we can use the Newton-Raphson method: given a 
function (hopefully smooth) f(R) which vanishes when 
Pt{R) does, an estimate of the nodal distance is: 



d(R) 



\W)\ 



This is good as long as the contribution of higher order 
derivatives of / is not large within a thermal wavelength. 

For the particle-in-a-box model, we took f(x) = 
sin(7ra/a) for the purposes of ascertaining the effect of 
estimating inaccurately the nodal distance. This gives 
an estimate of the nodal distance d = (a/ir)\ tan(7ra/a)|. 
The energy with this assumption is shown as the mid- 
dle dashed line: the error converges as r 1 5 . The action 
is worse because the precise distance of the walk from 
the node is not known. Nonetheless, the error is much 
superior to that using the primitive action. 

In practice, one should choose / differently in the low 
temperature and high temperature limits, as shown in 
figure (8). At low temperatures, as we have discussed, 
the nodes are smooth and relatively independent of R* 
and t. Then we choose f(R) = p?{R, R*',i)- On the 
other hand for small values of t the density matrix is very 
sharply peaked. The curvature of Gaussian density ma- 
trix always places the distance too close, so we consider 
an alternative procedure where one divides out the near- 
est Gaussian. The idea is to make / as linear as possible. 
Define f(R) = p T (R, R ;t)/ p(R,VR,;t) where V is the 
permutation having the largest term in the determinant: 
the permutation which minimizes X]i( r i* — r, p) 2 - To 
determine the best permutation, we use the Hungarian 
algorithm for solving the linear sum assignment problem 
(Burkard and Derigs, 1980). That takes 0[N 3 ] opera- 
tions. 

Although we have made lots of approximations in de- 
riving the nodal action, it is correct in the limit that 
d t — > 0. This is an essential property of any approxima- 
tion. It matters little what you do for particles far away 
from the nodal surface. At high temperature, the nodes 
are far away (it some sense that are as far away as they 
can be) so the upper bound is reasonable. At low tem- 
perature, the wave function and nodes are smooth and 
planar, see figure 4. There the Newton-Raphson proce- 
dure works well. 

Now consider the problem of how to calculate the in- 
ternal energy. The simplest way of calculating the energy 
is to take the time derivative of the action. This is called 
the thermodynamic estimator in Ceperley (1995). The 
direct fermion contribution to the energy of the nodal 
action is: 



dU 



N 



dr 



M I Xi f 



d[ln(didi-i)] 



dt 



}) 



(8.6) 



(8.5) 



where X{ = c?,c?,_i/r. This is only the direct contribu- 
tion. The restriction changes the distribution of paths, 
so all terms in the energy are affected by Fermi statistics. 
In the above formula, one can neglect the derivative of 
the nodal distance on "time" . This is appropriate at low 
temperature but could cause some errors at intermediate 
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temperatures. Note that in computing the energy, all the 
time slices contribute, not just the reference point, since 
all the points contribute to the action. 

Let us now look at how the restriction is calculated 
during the path integral Monte Carlo random walk. Con- 
sider the case where the trial function is a Slater determi- 
nant of single particle density matrices as in Eq. (7.4). 
Typically in the Monte Carlo random walk only a few 
(say p) particles are moved at once. This means only p 
columns of the Slater matrix are changed. Recalculation 
of the full determinant will take 0[iV 3 ] operations but an 
update of a p atoms is much more rapid as was shown 
by Ceperley, Chester and Kalos (1977). The matrix of 
g(r,t, rj»; t) = gijt and its transposed inverse c;jt, defined 
as the solution to: 

^giktCjkt = Sij (8.7) 

k 

are computed and saved throughout the run. Using c, 
we can determine very quickly if a proposed trial move 
has a positive trial density matrix, because c is also the 
cofactor matrix of g. 
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FIG. 8. The antisymmetric density matrix for reference point 
at 1 and the current point at 1/2. Two ways of estimating the 
distance to the nodes, both using the Newton-Raphson for- 
mula, are shown. In the first method the distance is inversely 
proportion to the logarithmically derivative of the density ma- 
trix, leading to an underestimate of the distance. In the sec- 
ond method, the estimate is the logarithmic derivative of the 
nearest neighbor is subtracted out. This leads to a distance 
estimate too large. 

Suppose a single particle with label i at time t is moved 



to a trial position: r, t — > r 7 . Then we compute the trial 
density matrix as: 

p T (R',R*;t) = p T (R, R f ;t)^2g(r' ,r kf ;t)c ikt . (8.8) 

k 

This takes 0[N] operations. If the trial density matrix 
is negative, the trial move is rejected. Otherwise, we 
go back and recompute the new distances and accept 
or reject on the resulting changed nodal action. The 
derivative needed for the distance in the nodal action 
is determined as: 

V'j In p T (R', R*;t) = J2 V</(r', i k .;t)c ikt (8.9) 

k 

Once a move is accepted, we must then update Ci kt . The 
Sherman-Morrison (1949) formula expresses how to cor- 
rect an inverse that is changed by a single column 

c 'jkt = c jkt + [Sji - bj]c ikt /bi (8.10) 

where bj = y^ k giktCjkt The update takes C[iV 2 ] opera- 
tions. After each update, we verify a few entries of Eq. 
(8.7). If the inverse is no longer accurate, the inverse 
is recomputed using Gaussian elimination with pivoting. 
This only happens every several thousand steps, but oc- 
curs more often as the number of particles increases or the 
temperature is lowered since overlapping orbitals make 
the Slater matrix ill-conditioned. Good update proce- 
dures for the more complicated backflow determinants 
have not been developed. It would be unfortunate if the 
more accurate trial density matrices were much slower to 
compute. 

B. Sampling restricted paths 

The other crucial technical aspect of path integrals is 
the question of how to sample the path integrals and 
how to do the permutational sum. While, in principle, 
one could use the Molecular Dynamics method for the 
path average, there seems to be no way of performing 
the permutational sum with MD, so in practice all calcu- 
lations where particle exchange is important have used a 
generalization of Metropolis Monte Carlo. In single slice 
sampling, a single atom in a single time-slice is moved and 
the new trial position is accepted or rejected based on the 
action of the new path. However, because the "beads" 
on the "polymer" are hooked together, uncorrelated mo- 
tion of individual beads is damped, so that the polymer 
as a whole moves very slowly. Multiple-slice sampling 
methods have been devised which speed up the motion 
through path space. These and the methods for mov- 
ing through permutation space are discussed at length in 
Ceperley (1995) for bosonic path integrals. 

Since the restricted path method is essentially built on 
top of the bosonic method, all of those techniques can 
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be used for fermion path integrals. Here we just remark 
on a few of the differences. First there is the question of 
moving non-reference point beads. We move them as if 
they were bosons, then once an acceptable bosonic path 
is found, we check the fermion action. First we check 
whether the new trial path has a positive action at all 
the new points. If it doesn't that path is immediately 
rejected. If it is positive, then we compute the estimated 
distance to the nodal surface, and the image action. This 
is used to finally accept or reject the path. 

The second difference is how to move through permuta- 
tion space. For bosonic problems, we try inserting cyclic 
permutations of two, three or four atoms at a time. For 
fermions we try only three atoms at a time since permu- 
tations for two or four atoms would certainly be rejected. 
A table of the free particle density matrix changes for all 
possible triplet permutations between particles with the 
same spin is constructed. This table is used to sample the 
permutation change. After the permutation is sampled, 
an appropriate path is sampled. 
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FIG. 9. The acceptance ratio as a function of temperature 
for liquid helium for a reference point relabeling move. In 
this move, another time slice is chosen at random to become 
the reference point. The move is accepted if the current path 
with the new reference point is node- avoiding. It rapidly ap- 
proaches zero at low temperatures. 

Finally, there is the question of how to move the ref- 
erence point. If the reference point is moved, the trial 
density matrix must be recomputed for all M — 1 links 
of the path. If any are found to be negative or too close 
to zero the reference point move is rejected. The action 
seen by the reference point is non-local. One might worry 



that this will trap the reference point, and the path as a 
whole. 

In a reference point relabeling move, one does not 
change the physical location of the paths but merely 
shifts the whole path in "time" . For example one could 
try to move the reference point from time t\ to time t?- 
The bosonic action is not affected by this shift so it does 
not have to be recomputed. One simply has to see if 
the path with the new reference point is node-avoiding. 
Rejections occur only because the nodes have a depen- 
dence on time. In figure 9 is shown the acceptance ra- 
tio of this reference relabeling moves in liquid 3 He using 
free-particle nodes. It drops quickly to zero once the sys- 
tem becomes a good Fermi liquid. The reference point 
strongly resists any uniform shift in the time direction. 

We have found that one can treat movement of the 
particles at the reference point just like the other other 
time-slices. One moves a few atoms to a new trial po- 
sition and accepts or rejects that move. In this case we 
have to check whether that causes the trial density matrix 
to become negative at some later time. However, if only 
a single atom is being moved at the reference point., the 
acceptance ratio does not go to zero at low temperatures. 

Trapping has been seen in a related method for the 
lattice Hubbard model. The constrained path algorithm 
is similar to restricted path integral Monte Carlo, except 
the paths are in determinant space, not in configuration 
space. The algorithms both have a non-local fermion ac- 
tion. Fahy and Hamann (1990) had considered a path 
integral version of the algorithm and had seen trapping 
and non-ergodic behavior. In a recent paper, Zhang, 
Carlson and Gubernatis (1995), showed that the zero- 
temperature diffusion Monte Carlo algorithm does not 
have this trapping. 

We believe that sampling techniques can significantly 
reduce the likelihood of trapping but clearly non-ergodic 
behavior is possible in principle at sufficiently low tem- 
peratures or for large numbers of particles. Zero temper- 
ature algorithms appear not to suffer from trapping, but 
the importance sampling which is introduced can bias the 
calculated result. More experience is needed to find out if 
trapping is a generic problem for finite temperature path 
integrals, and whether there are general solutions to this 
problem. 

IX. PERMUTATIONS, THE MOMENTUM DISTRIBUTION 
AND FERMI LIQUIDS 

Here we discuss further qualitative characterization 
of restricted paths. In our bosonic paper (Ceperley, 
1995) was discussed many relationships between paths 
and physical properties. Much less is known for re- 
stricted fermion path integrals but the nodal restriction 
and the extra spin degrees of freedom make them poten- 
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tially much richer. Since one can compute the fermion 
partition function with restricted paths and all thermo- 
dynamic properties can be obtained by differentiating the 
partition function, the paths should reflect any kind of 
ordering that the fermion system undergoes. 

Superfluidity in liquid 3 He was first detected by mea- 
suring the specific heat. It is interesting to speculate how 
this triplet pairing will appear with restricted fermions 
paths. As Feynman(1953) first pointed out, superfluidity 
in bosonic system is signaled in path integrals by the ap- 
pearance of permutation cycles of macroscopic size. This 
has a different interpretation for fermion paths: long per- 
mutations are related to the formation of a Fermi liquid, 
as we now describe. 

Let us consider how to calculate the momentum distri- 
bution with restricted paths. The momentum distribu- 
tion is the Fourier transform of the single particle density 
matrix. For a translationally invariant system: 



n(r) 



(k F r)- 



;[sin(kpr) — kprcos(kpr)] (9.6) 



8^ 3 J 



dre- ik - r n(r) 



(9.1) 



where the single particle density matrix is: 

n(r) = ^ J dRp{T X + r, r 2 , . . .,r n ,ri,r 2 , . ..,r N ;/3). 

(9.2) 

The momentum distribution is the Fourier transform of 
an off-diagonal element of the density matrix. The paths 
that we have been discussing up to this point, where all 
of the paths end at the start of another particle's path, 
cannot be used to calculate the momentum distribution. 
To get an observable in momentum space, we cannot do 
the simulation entirely in the position representation. All 
that is needed to get the momentum distribution is to 
allow one of the atoms to have free ends. 

However, now the fermion sign comes back in, as it 
must. First, consider the classical Maxwellian distribu- 
tion for which: 



? 'k \ ' ex P(-/ ?Ak2 ) 



n(r)=exp(-— ) 



(9.3) 
(9.4) 



Here the two ends of the cut polymer have an end-to-end 
vector r and typically will be within a thermal wavelength 
of each other. At high temperature, the polymers do not 
exchange and the shape of an individual atom's path is 
a pure Gaussian random walk. Hence, both the single 
particle density matrix and the momentum distribution 
are Gaussian; both smooth analytic functions. 

Now, consider the ideal fermi-gas momentum distribu- 
tion (for spinless fermions). 
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/(2tt 3 P ) for k < kp 
for k > kp 



where the Fermi wavevector for spinless fermion is related 
to the density by kp = (6tt 2 p) 1 ^ 3 . Note that the single 
particle density matrix is proportional to the spherical 
Bessel function ji(z), and slowly decays to zero at large 
r. It has zeroes at kpr = 4.493, 7.725, . . .. These zeroes 
mark the places where the even and odd permutations 
cancel out. Since n(r) is often negative, even with re- 
stricted paths we must have negative weights entering. 
The momentum distribution has a discontinuity at the 
Fermi wavevector lip. As a consequence the single par- 
ticle density matrix must decay at large distance as r -2 . 
We can get such long-range behavior only if there are 
macroscopic exchanges. Hence the existence of any kind 
of non-analytic behavior (we mean a discontinuity in 
or in any of its derivatives) implies that the restricted 
paths have important macroscopic permutation cycles. 

The calculation of the momentum distribution using 
restricted paths has not yet been performed. The Monte 
Carlo procedure would include a pair exchange between 
the cut end and another atom. Remember that for di- 
agonal simulations, only three-body permutations are al- 
lowed by the restriction. The pair moves generate the 
negative contributions to the single particle density ma- 
trix. 
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FIG. 10. The probability that an atom is on a permutation 
with n other atoms shown at three temperatures, the low- 
est curve with triangles is at 2K, with squares is at IK, with 
circles 0.5K. The horizontal line at 1/19 is the expected dis- 
tribution for a random permutation. 
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In figure 10 is shown the permutational cycle distri- 
bution for liquid 3 He, a typical Fermi liquid. It is seen 
that at temperatures below 0.5K the cycle distribution 
approaches a uniform distribution. Above this temper- 
ature, there is a substantial even-odd effect in the cycle 
distribution function. 

X. OTHER FERMION METHODS 

Let me briefly mention several other techniques to han- 
dle fermion path integrals in the continuum. 

A. Hall's method 

A very similar method to what we have described here 
is an approximate projector method(Hall, 1991, Hall and 
Price, 1991, and Hall, 1992). Hall argues that walks that 
are near an exact node of the density matrix end up not 
contributing to the value of the partition function (or any 
other property) but do contribute to the variance. Hall 
discards all paths for which \pT(R t , R*; ft)\ < e with the 
method becoming exact as e — > 0. Note that it is possible 
to get negative contributions; the walks can "hop" over 
the nodes. Hall argues that if £ >> r one gets rid of 
most of the negative signs, but the result is approximate. 
If c « r, hopping is easy, the method is exact but the 
variance is just as large as it was before the restriction. 
Hence, e and r are parameters that can be used to control 
independently the statistical error and the systematic er- 
ror. 

From the restricted-path point-of-view we can under- 
stand how Hall's method can work. If we eliminate paths 
close to a node, we eliminate equal numbers of positive 
and negative paths and this will not change the answer. 
But by formulating the restriction as a boundary condi- 
tion instead of as a condition on the value of the density 
matrix, we can make better nodal actions, and work at 
larger values of e and r without losing accuracy. 

Hall also symmetrizes the action with respect to which 
point is considered as the reference point. This will in- 
crease the computational time since M 2 checks of the 
trial density matrix have to be made versus M if you 
do not symmetrize. Our experience with liquid 3 He for 
temperatures below the degeneracy temperature, using 
free-particle trial density matrix, shows that only a sin- 
gle point at a time can be the reference point for a typical 
path. Symmetrizing with respect to the reference point 
only increases the computational work, without bring any 
increased efficiency. 

Hall's method may be useful in making a path integral 
version of the release-node method that is used in ground- 
state Quantum Monte Carlo to go beyond the fixed-node 
approximation. One would like to have a method to add 



in a few node-crossing walks to fix up the nodal surfaces 
locally, but in a controlled fashion so that the efficiency 
is no more than one order-of-magnitude lower than the 
restricted-path efficiency. 

B. Slice-wise antisymmetrization 

One is allowed to apply the antisymmetrization pro- 
jection operator on each time slice as suggested by Taka- 
hashi and Imada (1984). The action contains a determi- 
nant at each slice, not a restriction. The product of the 
signs of the determinants is used to weight the paths as 
discussed in Sec. (3). In contrast to the restricted path 
integral approach, one retains time slice symmetry. In 
one spatial dimension this solves the sign problem since 
one picks up a sign when two fermions hop over each 
other, but if one antisymmetrizes each step, the positive 
sign path always be greater than the negative path. Con- 
sider this procedure as applied to the hydrogen molecule 
in Sec. (5). The probability of a path will be the bosonic 
action times a factor: 

M 

PF <x 11 [1 - exp(-r k • r*_i/(Ar))] (10.1) 

where r = ri — r?- This prevents an exchange in a single 
step (actually it only rules outsteps with |rft-r^ + i| < At). 
But it does not rule out paths of type B although it may 
diminish their probability somewhat if the thermal wave- 
length is comparable to the bond-length. The antisym- 
metrization has a spatial range r ss (Ar) 1 / 2 so it is only 
a local modification of the action. At low temperature 
there will still be a cancellation between paths A and C 
against those of B resulting in an exponentially bad sign 
problem. 

Newman and Kuki(1992) have developed this tech- 
nique further by averaging also over other rotational sym- 
metries. They achieved much higher efficiencies for two 
electron problems than with the direct method. They 
still observed an exponentially decreasing efficiency, as /3 
increased, but the rate of growth was much diminished. 
With increasing numbers of particles the effect of this 
symmetrization is less. If the boson system can undergo 
a lambda transition, the fermion efficiency is guaranteed 
to be exponentially small. The conclusion is that single- 
slice antisymmetrization will not help much for interest- 
ing fermion systems. Also, it does not seem easy to com- 
bine it with the restricted path integral approach. One 
has to choose between antisymmetrizing at every slice, 
or restricting the paths to be node-avoiding. 
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C. Ignoring the sign 

Suppose one does as much slicewise antisymmetriza- 
tion as possible but ignore the sign of the negative walks. 
This was suggested by Imada (1984), giving the following 
ingenious statistical argument: 

In principle, even one sample [of an in- 
tegrand] can give a good convergence to the 
exact value if the calculation is performed 
in the thermal equilibrium state of a suffi- 
ciently large system. It should be noted that 
the thermal equilibrium state is realized in 
a sample generated by the importance sam- 
pling method, in which the configurations are 
generated in proportion to the contribution to 
the partition function. If we regard the gener- 
ated samples of a finite system as a sequence 
in time evolution, our assertion seems to be 
related to the ergodic hypothesis: A long time 
average is equal to the phase average. If we 
calculate the average from only one sample 
Eq. (3.2) is replaced by 

(A) £3 X\ f% Xi ^ A3 . . . f% Xq. (10.2) 

If the system size is infinite, X{ should not 
depend on i. When the system size is large 
but finite in the actual simulation, better con- 
vergence is expected from the average: 

G 

{X)^J2x i /G. (10.3) 

i = l 

The finite size effect in Eq. (10.3) must de- 
crease with increasing system size. Equation 
(10.3) is a surprisingly simple result. This 
can be calculated without concern regarding 
the sign of the weight function. As similar to 
the ergodicity in classical systems, the valid- 
ity of Eq. (10.3) in sufficiently large systems 
seems to be obvious from a physical point of 
view, although it has not been mathemati- 
cally proven yet. 

He concludes that ignoring the sign is an exact procedure 
for a sufficiently large system. 

Such a procedure cannot be correct. As we have ex- 
plained above, when one ignores the sign, one is using an 
effective bosonic action. If the system remains liquid as 
it is cooled, it must Bose condense. One will get a com- 
pletely different behavior in a fermion system. The above 
reasoning is flawed because classical statistical mechanics 
applies to a probability distribution, not to an integrand 
which can be positive and negative. [See also Sorrella et 
al (1988,1989) for a related method for lattice models and 
the analysis of the Sorrella method by Loh et al (1990).] 



D. Cancellation 

Cancellation schemes have been developed for zero 
temperature quantum Monte Carlo, see for example An- 
derson (1995). By cancellation, we mean that one starts 
with a population of positive and negative random walks. 
If a positive and negative random walk approach each 
other, they annihilate. Alternatively, they could repel 
each other. This formally solves the sign problem. The 
remaining sticky point is to determine how many walks 
are needed to create a stable interface between the pos- 
itive and negative walks. Estimates show that the min- 
imum stable population size will grow exponentially in 
the number of fermions. Hence, this has not yet led to a 
viable many-fermion algorithm. Such methods have not 
yet been developed for finite-temperature path integrals. 

XI. CURRENT APPLICATIONS OF FERMION PATH 
INTEGRALS 

I will not make a detailed comparison between the 
path integral approach for quantum simulations with the 
other methods such as Variational Monte Carlo, projec- 
tor Monte Carlo (Greens function or Diffusion) or ab ini- 
tio molecular dynamics using forces from LDA. There 
is a comparison between boson path integrals and zero 
temperature quantum Monte Carlo methods in Ceper- 
ley (1995). Path Integrals are unique in being able to 
treat strongly correlated systems, and mixed quantum- 
classical systems, including the effects of temperature, 
and having other, more technical advantages. 

I also think that they are unique in being amenable 
to a "black-box" approach since no trial function is used, 
only a trial restriction. The best proof of this assertion is 
that applications of RPIMC to a variety of systems are al- 
ready beginning to appear. Some of the applications are 
considerably more complicated that have been treated 
with zero-temperature quantum Monte Carlo methods. 

The first application of RPIMC was to liquid 3 He by 
Ceperley (1992) using the semi-empirical Aziz potential. 
Using free-particles nodes, these calculations showed rea- 
sonable agreement with experiment and also showed the 
importance of using time-dependent nodes and permuta- 
tions. See figure 11. There has also been a recent study 
of the interface between superfluid 4 He and liquid 3 He 
by Boninsegni and Ceperley (1995) who were able to de- 
termine the width of the interface and estimate directly 
solubilities and effective masses in the 3 He-rich and 4 He 
rich phase. 

In our group, there has been a major effort to calcu- 
late properties of an electron-proton hydrogenic plasma 
by Pierleoni et al. (1994) and Magro et al. (1994). In 
these calculations, from 32 to 64 electrons and an equal 
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number of protons were put into a periodic box, interact- 
ing only with the Coulomb potential. Both the electrons 
and protons were fully quantum particles. A time step 
of roughly 10 6 K was found adequate to give accurate re- 
sults. We were thus able to reach temperatures of 4000K 
using 250 time slices. Using reasonable length runs on 
CRAY-YMP computers and fast workstations we found 
good agreement with other theoretical approaches at high 
density and at high temperature in the plasma phase. 




T (K) 



FIG. 11. The energy of liquid 3 He vs temperature at a density 
of 0.16355 A -3 (near zero pressure). The line represents the 
experimental data. The open symbols are the RPIMC calcula- 
tions using ground state nodes; the filled symbols are RPIMC 
calculations with the free-particle density matrix nodes. The 
circles are with permutations, the triangles are with the iden- 
tity permutation. The zero temperature result was obtained 
with GFMC with free-particle nodes. 

At temperatures below 10,000K we observed the spon- 
taneous formation of molecules, as evidenced by a strong 
peak in the proton-proton correlation function at a dis- 
tance of lA. At slightly lower temperature of 7,000K we 
observed a very rapid build up in the molecular den- 
sity, enough to suggest that there is a first-order phase 
transition from the plasma state into a molecular liquid. 
The existence of the so called "plasma phase transition" 
has been guessed at for a few years. It can be crucial 
for understanding the interior of the giant planets or, 
more generally, the phase diagram of hydrogen at high 
temperatures and pressures. The ability to study such 
phase transitions at finite temperature shows the power 
of the new method. Also striking is the lack of theoreti- 
cal input to the calculation. Equivalent Variational and 



Greens Function Monte Carlo calculations at zero tem- 
perature (Ceperley and Alder, 1987, Natoli et al, 1993, 
1995) required careful construction and optimization of 
trial wave functions. 

These various simulations demonstrate the power of 
the method. Once the RPIMC method is programmed, 
one can rather directly do highly accurate simulations 
of experimentally relevant systems without tedious con- 
struction of basis functions and trial functions; the type 
of systems which can be treated are much more complex 
than can be handled with other techniques. 

XII. FUTURE PROSPECTS 

We think that the Path Integral Monte Carlo method 
is a very powerful method capable of doing very accu- 
rate simulations of both boson and fermion systems. Its 
methodology is much less well developed than for clas- 
sical systems and there are many challenges in applying 
this method to real systems. Among them are: 

1. Clearly much work needs to be done in figuring out 
what we should use as nodes since the restriction is 
the only uncontrolled approximation. Free-particle 
nodes are too primitive and possibly pathological 
in that they do not have any coupling from un- 
like spins. Analytical and numerical methods are 
needed to go from the free particle nodes at high 
temperature to the mean-field-like nodes that are 
useful at low temperature for atoms and molecules. 
We are currently working on this. We need to clar- 
ify the situation with regard to free-energy bounds 
so as to internally (within the RPIMC method) de- 
termine optimal nodes. 

2. We need ways to get to lower temperatures. One of 
the advantages of the RPI approach is that systems 
with both quantum and classical degrees of freedom 
are treated more easily than in other approaches. 
But to exploit this possibility, we need to simulate 
electrons at room temperature. So far, our lowest 
temperature with many-body hydrogen including 
both electron and proton paths has been 3,000K. 
There are several ways of going to lower tempera- 
ture, for example one could use zero- temperature 
nodal surfaces for the electrons. But much experi- 
ence is needed to verify these approximations. 

3. One of the big advances in recent years in zero- 
temperature quantum Monte Carlo methods has 
been in the development of pseudo-potential meth- 
ods. This has allowed one to treat systems with 
carbon, nitrogen, silicon atoms as easily as hy- 
drogen and those with transition metal atoms (e. 
g. iron and copper) with an order of magni- 
tude more effort. Two basic approaches have been 



21 



developed. In the pseudo-Hamiltonian approach 
(Bachelet, Ceperley and Chiocchetti, 1989), one al- 
ters the electron mass inside the core so one can 
use directly all the quantum Monte Carlo methods. 
With this approach, one cannot treat all types of 
atoms sufficiently accurately. In the more general, 
non-local pseudo-potential method one has a dif- 
ferent potential for each angular momentum state 
of the electron. The angular momentum depen- 
dence causes problems with the nodal restriction 
since paths are not continuous. What has been 
done (Mitas et al 1991) is to project the non-locality 
onto an accurate trial function and thereby make 
it local. One only requires accurate trial functions 
in the core region. Perhaps the use of atomic wave 
functions will be sufficient to make this work in 
path integrals. We do not want to have to bring 
many-body trial functions back in. 

4. Path Integral Monte Carlo methods are computa- 
tionally expensive. This has limited the number 
of applications. The easy answer to this objection 
is that computers are getting faster and cheaper. 
Nonetheless, improving the efficiency of the meth- 
ods will allow us to solve much more difficult or re- 
alistic problems. More effort needs to go into ways 
of speeding up the dynamics of moving the paths, 
particularly at low temperature. We also need bet- 
ter nodal actions so the time step can be made 
larger. Finally, we need to write C[N] codes so we 
can do much larger systems. (Note that RPIMC 
is naturally C[N] since the single-electron density 
matrix in Eq. (7.5) is localized in the range of the 
thermal wavelength.) 

5. The computation of dynamical properties is a chal- 
lenge. As Gubernatis has described there is much 
progress with maximum entropy methods. RPIMC 
has gotten rid of the sign problem in calculating the 
partition function but what is needed for dynam- 
ics are off-diagonal density matrix elements. We 
have gotten rid of most but not all of the nega- 
tive permutations on the diagonal (Boninsegni and 
Ceperley 1994). Because of these signs, it seems 
unlikely that RPIMC will be able to calculate the 
highly accurate imaginary-time correlation func- 
tions needed for maximum entropy reconstruction 
of the real-time dynamics. As a consequence, prop- 
erties such as conductivity, and spectral functions, 
seem very difficult to compute with these methods. 
Perhaps one needs to formulate dynamical proper- 
ties in terms of thermodynamic properties as done 
with the superfluid density of bosonic systems. 

6. Much needs to be done in understanding how to in- 
terpret RPI. For example, what is the description of 



superfluidity, the Mott transition, molecular bond- 
ing and bound versus continuum states? Girardeau 
(1990) has suggested that eigenvalues and eigen- 
functions of the one- and two-body density matri- 
ces can be used to identify the number of molecules. 
Such eigenfunctions (the natural orbitals) are com- 
putable with path integrals as we described above 
with the momentum distribution. This shows how 
to make the connection between properties of the 
paths and rigorous many-body theory. There are 
a few hints that many of the relationships we have 
discovered for bosonic systems have equally good 
interpretations in terms of RPI. 

7. The generalization of the method to magnetic fields 
is possible, see Ortiz et al (1993). It will also be in- 
teresting to work with other ensembles, for example 
with fixed angular momentum to see vortices. 

The above is a list of challenges for those who wish 
to develop further fermion path integrals. Hopefully, we 
will have some answers by the next simulation meeting 
at Lake Como. 
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